home *** CD-ROM | disk | FTP | other *** search
/ Gold Medal Software 3 / Gold Medal Software - Volume 3 (Gold Medal) (1994).iso / stats / lsqrft15.arj / MRQFIT2.C < prev    next >
Text File  |  1994-01-01  |  8KB  |  303 lines
  1. #include <math.h>
  2. #include <stdio.h>
  3. #include <stdlib.h>
  4. #include "fit.h"
  5.  
  6. /* calculates alpha, beta, and chisq */
  7. /* see Numerical Recipes in C for definition of these */
  8. /* the Levenberg Marquardt algorithm requires us to slove */
  9. /* the equation alpha*da = beta (where alpha is a matrix */
  10. /* and da and beta are vectors) for da */
  11. /* this function does this by gauss-jordan */
  12. /* elimination */
  13.  
  14. /* declarations for data declared externally */
  15. extern int grflag;
  16. extern int gnuopen;
  17. extern FILE *pstream;
  18. extern double *xmin, *xmax;
  19. extern int wiflag;
  20. extern int veflag;
  21. extern int debug;
  22. extern char *GNUPLOT;
  23. extern int MYPLOT;
  24. extern int doflag;
  25.  
  26. int mrqfit(data,order,num_indep,ndata,itmax,a,ma,lista,mfit,
  27.                 covar,chisq,func,filename,comment)
  28. double **data;   
  29. struct data_order order;
  30. int num_indep;
  31. int ndata;
  32. int itmax;
  33. double *a;
  34. int ma, *lista, mfit;
  35. double **covar, *chisq;
  36. int (*func)();
  37. char *filename;
  38. char *comment;
  39. {
  40. int i, j, k;  /* indices for loops */
  41. int failed;
  42. double **alpha;
  43. /* alamda determines how much we change the fitting parameters */
  44. double alamda;
  45. char cmd[20]="";
  46. double *beta;
  47. double *atry; /* parameters to try and see if we reduce chisqr */
  48. double **alpha_try;
  49. double ochisq; /* best chisqr so far */
  50. double *da;
  51. FILE *stream;
  52.  
  53. /* allocate space for arrays */
  54. atry = dvector(ma);
  55. alpha_try=dmatrix(ma,ma);
  56. beta = dvector(ma);
  57. da = dvector(ma);
  58. alpha=dmatrix(ma,ma);
  59.  
  60. stream=NULL;
  61.  
  62. /* calculate alpha, beta, and chisq for current value of parameters */
  63. if(debug) printf("in mrqfit, calling alpha_beta_chisqr()\n");
  64. failed=alpha_beta_chisq(data,order,num_indep,ndata,a,ma,lista,mfit,
  65.                             alpha,beta,chisq,func);
  66. if(debug) printf("in mrqfit, alpha_beta_chisqr() returned %d\n", failed);
  67. if(failed){
  68.     free_dmatrix(alpha,ma,ma);
  69.     free(atry);
  70.     free_dmatrix(alpha_try,ma,ma);
  71.     free(beta);
  72.     free(da);
  73.     return 0;
  74. }
  75.  
  76. /* start alamda small */
  77. alamda = 1e-3;
  78. i = 0;
  79.  
  80. /* loop until user tells us to quit */
  81. while( strcmp(cmd,"q") != 0){
  82.  
  83. /* loop specified number of iterations */
  84. while( i <= itmax ){
  85.     ochisq = *chisq;
  86.     for(j=0;j<mfit;j++){
  87.         for(k=0;k<mfit;k++){
  88.             alpha_try[j][k] = alpha[j][k];
  89.         }
  90.         alpha_try[j][j] = alpha[j][j] + alamda;
  91.     }
  92.     if(debug)printf("in mrqfit, calling solve_for_da()\n");
  93.     solve_for_da(alpha_try, covar, beta, da, mfit);
  94.     if(debug)printf("in mrqfit, returned from solve_for_da()\n");
  95.  
  96.     for(j=0;j<ma;j++) atry[j] = a[j];
  97.     for(j = 0; j < mfit; j++) atry[lista[j]] += da[j];
  98.  
  99. if(veflag==2){
  100.         printf("\nalpha: ");
  101.         for(j=0;j<mfit;j++){
  102.             printf("\n %d   ",lista[j]);
  103.             for(k=0;k<mfit;k++) printf(" %g ",alpha_try[j][k]);
  104.         }
  105.         printf("\ncovar: ");
  106.         for(j=0;j<mfit;j++){
  107.             printf("\n %d   ",lista[j]);
  108.             for(k=0;k<mfit;k++) printf("\n covar[%d][%d]= %g ",lista[j], lista[k], covar[j][k]);
  109.         }
  110.         printf("\natry:");
  111.         for(j=0;j<mfit;j++) printf("\n atry%d = %g",lista[j], atry[lista[j]]);
  112.         printf("\n");
  113.     }
  114.  
  115.     if(debug)printf("in mrqfit, calling alpha_beta_chisqr()\n");
  116.     failed=alpha_beta_chisq(data,order,num_indep,ndata,atry,ma,
  117.                             lista,mfit,alpha,beta,chisq,func);
  118.     if(debug)printf("in mrqfit, alpha_beta_chisqr() returned %d\n", failed);
  119.  
  120. if(failed){
  121.     free_dmatrix(alpha,ma,ma);
  122.     free(atry);
  123.     free_dmatrix(alpha_try,ma,ma);
  124.     free(beta);
  125.     free(da);
  126.     return 0;
  127. }
  128.  
  129.     if(*chisq >= ochisq){
  130.         alamda *= 10;
  131.         if(alamda > 1e15) alamda /= 3e14;
  132.         *chisq = ochisq;
  133.     }
  134.     else{
  135.         if(debug){
  136.             if(stream == NULL){
  137.                 printf("opening lamda.sts\n");
  138.                 stream = fopen("lamda.sts","at");
  139.             }
  140.             if(stream != NULL){
  141.                 fprintf(stream,"%g\n",alamda);
  142.                 printf("alamda: %g\n",alamda);
  143.             }
  144.         }
  145.         alamda *= 0.1;
  146.         if(alamda < 1e-15) alamda *= 3e14;
  147.         for(j=0;j < mfit; j++)a[lista[j]] = atry[lista[j]];
  148.     }
  149.  
  150.     if(veflag){
  151.         for(j=0; j < ma; j++)printf("\n a%d = %g", j, a[j]);
  152.         printf("\n chisqr = %g", *chisq);
  153.         printf("\n alamda = %g", alamda);
  154.         printf("\n");
  155.     }
  156.     else printf("iteration: %d chisqr: %g\n", i, *chisq);
  157.     i++;
  158.  
  159. #ifdef OS2
  160.     if( grflag && MYPLOT && (num_indep == 1) && veflag){
  161.         if(debug)printf("in mrqfit, calling myplot()\n");
  162.         failed=myplot(func,data, order, num_indep, ndata, filename, comment, a, ma);
  163.         if(debug)printf("in mrqfit, myplot() returned %d\n", failed);
  164.     }
  165. #endif
  166. }
  167.  
  168. if(grflag && ( !MYPLOT || !veflag || (num_indep != 1) )){
  169.         if(debug)printf("in mrqfit, calling plot()\n");
  170.         failed=plot(func,data, order, num_indep, ndata, filename, comment, a, ma);
  171.         if(debug)printf("in mrqfit, plot() returned %d\n", failed);
  172. }
  173.  
  174. if(doflag) strcpy(cmd, "q");
  175. else{
  176.     printf("q quit, g graphs, n does not graph, anything else continues ");
  177.     fflush(stdout);
  178.     gets(cmd);
  179. }
  180. i = 0;
  181. if(strcmp(cmd,"g") == 0) grflag = 1;
  182. if(strcmp(cmd,"n") == 0) grflag = 0;
  183.  
  184. }
  185. free_dmatrix(alpha,ma,ma);
  186. free(atry);
  187. free_dmatrix(alpha_try,ma,ma);
  188. free(beta);
  189. free(da);
  190. if(stream != NULL) fclose(stream);
  191.  
  192. return 0;
  193. }
  194.  
  195.  
  196. int alpha_beta_chisq(data, order,num_indep,ndata,a,ma,
  197.                                 lista,mfit,alpha,beta,chisq,funcs)
  198. double **data;
  199. struct data_order order;
  200. int num_indep;
  201. int ndata;
  202. double *a;
  203. int ma, *lista,mfit;
  204. double **alpha, *beta, *chisq;
  205. int (*funcs)(double *,double *,double *, double *,int, int, 
  206.                     int *, int *, double *, double);
  207. {
  208. int i,j,k;
  209. int failed;
  210. double *x, *y, *sig,*xsig;
  211. double ymod, *dyda, sig2i, wt,dy;
  212. int *fita, *dydx_flag;
  213. double *dydx;
  214. double sigi;
  215.  
  216. if(debug)printf("in abc, allocating data\n");
  217. dyda=dvector(ma);
  218. fita = ivector(ma);
  219. x = dvector(num_indep);
  220. dydx_flag = ivector(num_indep);
  221. dydx = dvector(num_indep);
  222.  
  223. if(debug)printf("in abc, setting pointers for y and sig\n");
  224. y = data[order.y];
  225. sig = data[order.sig];
  226.  
  227. if(debug)printf("in abc, setting pointers for xsig\n");
  228. for(i = 0; i < num_indep; i++){
  229.     if(order.xsig[i] >= 0){
  230.         xsig = data[order.xsig[i]];
  231.         dydx_flag[i] = 1;
  232.     }
  233.     else dydx_flag[i] = 0;
  234. }
  235. if(debug)printf("in abc, done setting pointers for xsig\n");
  236.  
  237. for(i = 0; i < ma; i++){
  238.     fita[i] = 0;
  239.     for(j = 0; j < mfit; j++) if(lista[j] == i) fita[i] = 1;
  240. }
  241.  
  242. *chisq = 0;
  243. for(j=0; j < mfit;j++){
  244.     beta[j] = 0;
  245.     for(k=0;k<mfit;k++) alpha[j][k] = 0;
  246. }
  247.  
  248. if(debug)printf("in abc, setting pointers for x\n");
  249. for(i=0;i<ndata;i++){
  250.     for(j = 0; j < num_indep; j++){
  251.         x[j] = data[order.x[j]][i];
  252.     }
  253. if(debug > 1)printf("in abc, done setting pointers for x\n");
  254.  
  255.     if(wiflag == 0 || window(x, num_indep)){
  256.         if(debug > 1)printf("in abc, about to call funcs\n");
  257.         failed = (*funcs)(x,a,&ymod,dyda,ma,1, fita, dydx_flag, dydx, y[i]);
  258.         if(debug > 1)printf("in abc, returned from funcs\n");
  259.         if(failed){
  260.         free(dyda);
  261.         free(fita);
  262.         free(x);
  263.         free(dydx_flag);
  264.         free(dydx);
  265.         return 1;
  266.         }
  267.         data[order.yfit][i] = ymod;
  268.         sig2i = (sig[i]*sig[i]);
  269.         for(j = 0; j < num_indep; j++)
  270.             if(order.xsig[j] >= 0) sig2i += (dydx[j])*(dydx[j])*xsig[i]*xsig[i];
  271.         if(debug==2) printf("i: %d x[0]: %g ymod: %g y: %g sig2i: %g\n", i, x[0], ymod, y[i], sig2i);
  272.         dy = y[i] - ymod;
  273.         for(j=0;j<mfit;j++){
  274.             wt=dyda[lista[j]]/sig2i;
  275.             for(k=0;k<=j;k++)alpha[j][k] += wt*dyda[lista[k]];
  276.             beta[j] += dy*wt;
  277.         }
  278.         *chisq += dy*dy/sig2i;
  279.     }
  280. }
  281.  
  282. for(j=1;j<mfit;j++)
  283.     for(k=0;k<=j-1;k++) alpha[k][j]=alpha[j][k];
  284.  
  285. free(dyda);
  286. free(fita);
  287. free(x);
  288. free(dydx_flag);
  289. free(dydx);
  290. return 0;
  291. }
  292.  
  293. int window(double *x, int num_indep){
  294.  
  295. int i;
  296.  
  297. for(i = 0; i < num_indep; i++)
  298.     if(x[i] < xmin[i] || x[i] > xmax[i]) return 0;
  299.  
  300. return 1;
  301. }
  302.  
  303.